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ABSTRACT 

Many spacecraft attitude determination methods use exactly two vector measurements. The two vectors are typically the unit 
vector to the Sun and the Earth’s magnetic field vector for coarse “sun-mag” attitude determination or unit vectors to two 
stars tracked by two star trackers for fine attitude determination. TRIAD, the earliest published algorithm for determining 
spacecraft attitude from two vector measurements, has been widely used in both ground-based and onboard attitude 
determination. Later attitude determination methods have been based on Wahba’s optimality criterion for n arbitrarily 
weighted observations. The solution of Wahba’s problem is somewhat difficult in the general case, but there is a simple 
closed-form solution in the two-observation case. This solution reduces to the TRIAD solution for certain choices of 
measurement weights. This paper presents and compares these algorithms as well as sub-optimal algorithms proposed by 
Bar-Itzhack, Hannan, and Reynolds. Some new results will be presented, but the paper is primarily a review and tutorial. 

INTRODUCTION 

Suppose that we have measured two unit vectors and b 2 in the spacecraft body frame. These can be the unit vectors to an 
observed object like a star or the Sun, or some ambient vector field such as the Earth’s magnetic field. We consider only unit 
vectors because the length of the vector has no information relevant to attitude determination. Each of these unit vectors thus 
contains two independent scalar pieces of attitude information. The spacecraft attitude is represented by a 3x3 orthogonal 
matrix A, Le. A T A = / , the 3x3 identity matrix. The attitude matrix must also be proper, i.e., it must have unit determinant, so 
it is an element of the three- parameter group SO(3). Euler’s Theorem states that the most general motion of a rigid body with 
one fixed point is a rotation about some axis. This shows explicitly that SO(3) is a three-parameter group, since the three 
parameters can be taken as the rotation angle and two parameters specifying a unit vector along the rotation axis. Thus two 
unit vector measurements determine the attitude matrix, in general; in fact they overdetermine it. 

It is also necessary to know the components of the two measured vectors r [ and r 2 in some reference frame. The reference 
frame is usually taken to be an inertial frame, but this is not necessary. One can use a rotating frame such as the frame 
referenced to the orbit normal vector and the local vertical. The attitude matrix to be determined is the matrix that rotates 
vectors from the reference frame to the spacecraft body frame. Thus we would like to find an attitude matrix such that 


Ar, = b, 

(la) 

and 


Ar 2 = b 2 . 

(lb) 

This is not possible in general, however, for equation (1) implies that 


b, b 2 =(/4r,)- (Ar 2 ) = r, r A r Ar 2 = r, r r 2 = r, • r 2 . 

(2) 


This equality is true for error-free measurements, but is not generally true in the presence of measurement errors. It will be 
seen in the following that all reasonable two-vector attitude determination schemes give the same estimate when equation (2) 
is valid. 


It is clear from simple counting arguments that the two independent scalar pieces of information contained in a single vector 
measurement cannot determine the attitude uniquely. More concretely, if the attitude matrix A obeys equation (la), then so 
does the matrix /?(b 1 ,0 A )A/?(r,,0 r ), for any <t> b and <p r , where R(e,<p) denotes a rotation by angle 0 about the axis e. This line 
of argument also makes it clear that the attitude matrix is not uniquely determined if either the pair ^ and b : or the pair r x 
and r 2 are parallel or antiparallel. 

The earliest published algorithm for determining spacecraft attitude from two vector measurements was the TRIAD 
algorithm 1 This algorithm has been widely used in both ground-based and onboard 3 attitude determination. The two vectors 
are typically the unit vector to the Sun and the Earth’s magnetic field vector for coarse “sun-mag” attitude determination or 



unit vectors to two stars tracked by two star trackers for fine attitude determination. Recent developments in star tracker 
technology have produced star trackers that can track 5, 6, or even 50 stars at a time. For attitude determination using more 
than two vectors, optimal estimators based on a loss function introduced by Wahba are appropriate 4 . However, Bronzenac 
and Bender have shown that the n vectors from a small-field-of-view star tracker can be replaced by an average vector 
without significant loss of precision 5 . With this approximation, the two star tracker case, even with multiple stars tracked in 
each star tracker, can be treated as a two-vector-measurement problem. 

With this motivation, we survey solutions to the two- vector measurement problem, beginning with TRIAD. We then consider 
the optimal solution of Wahba’s problem. After this, we look at sub-optimal algorithms have been proposed by Bar-Itzhack 
and Harman 6 and by Reynolds 7 ' 8 . We compare the various algorithms for both accuracy and computational effort, and finally 
present conclusions. 

TRIAD 

The TRIAD algorithm, introduced by Black in 1964 1,2 , is based on the following idea. If we have an orthogonal right-handed 
triad of vectors {v^ v 2 , v 3 } in the reference frame, and a corresponding triad {w ly w 2 , w 3 } in the spacecraft body frame, the 
the attitude matrix 

A =[w l :w 2 :w 3 ][v 1 :v 2 :v 3 ] 7 = w,v 7 + w 2 v[ + w 3 v[ (3) 


will transform the v, to the by 


A\ i = w, , i=l, 2, 3. 


(4) 


The TRIAD algorithm forms the tnad ( v,, v 2 , v 3 } from rj and r 2 , and the triad {w^ w 2 , w 3 } from b { and b 2 . Incidentally, 
TRIAD can be considered either as the word "‘triad 1 ’ or as an acronym for “TRIaxial Attitude Determination.” The triads can 
be formed in three convenient ways. First, it is useful to define the normalized cross products 

r 3 =(r, xr 2 )/|r, xr 2 | (5a) 

and 

b 3 s (b, xb 2 )/|b, x b 2 |. (5b) 

We note that r 3 or b 3 is undefined if the reference vectors or the observed vectors, respectively, are parallel or antiparallel. 

This is the case noted above in which there is insufficient information to determine the attitude uniquely. If this is not the 
case, two of the TRIAD attitude estimates are 


and 


A n = b,r 7 + b 3 r 3 7 + (bj x b 3 )(r x x r 3 ) 7 
A r2 =b,r 7 + b 3 r 3 7 +(b 2 xb 3 )(r 2 xr 3 ) 7 . 


and 


A T1 r x s b 2 (r, ■ r 2 ) + (b 2 x b 3 )[(r 2 x r 3 ) • r, ] = (r, • r 2 )b 2 + [b, - (b, • b 2 )b 2 ]|r, x r,|/|b, x b, . 



(6) 


(7) 

= b 2 , but 


| r i x r 2 |/|b, x b 2 | 

(8) 

|r, xr 2 /|b, xb 2 |. 

(9) 


Thus the estimate A n emphasizes the first measurement and A n emphasizes the second. It’s not difficult to see, though, that 
both A n and A n satisfy equations (la) and (lb) if b, • b 2 = r, ■ r 2 . 

The third form of TRIAD treats the two measurements symmetrically. We define the unit vectors 

r + s (r 2 4 - r, )/|r 2 + r,| = (r 2 + r, )/^2(l + r, • r 2 ) , 


r - = ( r : ' »*, )/|r, - r,| = (r 2 - r,)/^2(l - r, • r 2 ) , 


( 10 ) 

(ID 


and b + and b_ similarly. It is easy to see that r + is perpendicular to r_, b + is perpendicular to b_, and also that r 3 = r + x r_ and 
b 3 = b + x b_. Thus { r + , r_, r 3 } and { b + , b_, b 3 } are orthogonal triads, and the third TRIAD estimate is given by 

A 73 = b + r 7 + b_r 7 +(b,. x b J(r + x r_) 7 . 


( 12 ) 



This estimate treats the two observations symmetrically, and gives A T? r + = b + and A n r_ = b_, but 


A n r, = b t (r t T,) + b_(r_ •!•,) = - 


and 


A T 3 r z = M r + -r 2 ) + b_(r_ r 2 ) = 


1 + r s ■ r, 

lib, b. 


1 + r, • r 2 

1 + b| b 2 


(bj + b 2 ) + 


(bj + b 2 ) — 


1 r ' r - (b, -b,) 
I - b. • b, 1 - 


1 - r, • r, 
l“b, b 2 


(b,-b 2 ) 


Again, it’s not difficult to see that A n satisfies equations (la) and (lb) if b, b 2 = r { *r 2 . 


(13) 


(14) 


All three TRIAD estimates satisfy A n r 3 = b 3 , for i = 1, 2, 3. From this and the above observations, it is clear that A n , A n , and 
A-n give identical estimates if equation (2) is valid, since they provide the same mapping of a basis {r^ r 2 , r 3 } in the reference 
frame to a basis { bp b 2 , b 3 } in the spacecraft body frame. 

THE OPTIMAL SOLUTION 


In 1965, Grace Wahba, then a graduate student at Stanford University on a summer job with IBM, proposed the following 
problem 4 : Find the orthogonal matrix A with determinant +1 that minimizes the loss function 

< 15 ) 

where {bj is a set of n unit vectors measured in a spacecraft’s body frame, {r,} are the corresponding unit vectors in a 
reference frame, and {a,} are non-negative weights. We can rewrite equation (15), using the invariance of the trace under 
cyclic permutations, as 

= + l r i| 2 )-Z, a / b M r . =(X, a .)- trace ( /iBr )* ( I6 > 

where 

It is obvious that the attitude matrix that minimizes the loss function is the proper orthogonal matrix that maximizes 
trace(A5 r ). Almost all solutions of Wahba’s problem are based on this observation. The original solutions solved for the 
attitude matrix A directly, but most practical applications have been based on Davenport’s ^-method 2,9 , which solves for the 
attitude quaternion 10,11 . Shuster’s QUEST algorithm, in particular, has been widely used 12 . Shuster showed a simplification in 
the two-observation Wahba problem, but the first explicit closed-form solution was presented in reference 13. 

We begin by noting that the matrix B has the singular value decomposition 14 ' 16 

B = USV t , (18) 

where U and V are orthogonal matrices, and S is diagonal; 

S = diag(j p 5 2 ,s 3 ) , (19) 

with 

5 , > > ^ 3 > 0 . ( 20 ) 

In the two-observation case, it is clear from equation (17) that B has rank at most 2, and therefore 

det B = 5[J 2 5 3 = 0 . (21) 


Equations (20) and (21) show that 


j 3 = 0 


in the two-observation case. We shall take advantage of some resulting simplifications in this case. The general 
^-observation case is treated in references 13 and 14. 


( 22 ) 


Since s 3 = 0, we are free to choose the sign of the last column of U and of V so that both of these matrices have positive 
determinants. We shall assume that this is the case. Now 


trace(AB r ) = trace(.4V'S(/ r ) = trace(W'S) = s t W n + s 2 W n , 


(23) 



(24) 


where we have again used the invariance of the trace under cyclic permutations, and 

W = U t AV. 

Now using the Euler axis/angle parameterization for W = R(e,<p) gives' 011 

trace(AB r ) = j,[cos<p + e 2 (i - costp)] + s,[cos<p + e\(\ -cos^p)] = + J 2 e 2 + cos<p[-y,(l - e 2 ) + s 2 (l - e;)]. (25) 

This is clearly maximized for cos0 = 1, which means that W — l. Thus the optimal attitude is given by 

A op ^UV T . (26) 

Equation (25) shows that the minimum of traced# 7 ") is unique unless s, = 0. The vanishing of s 2 is the sign in the optimal 
algorithm that the observations are not sufficient to determine the attitude. We shall see below that this is related to the 
parallelism of the reference frame or body frame vectors. 

The singular value decomposition is rather expensive computationally, so we look for a simpler way to compute A opr We note 
that the classical adjoint, or adjugate, (the transposed matrix of cofactors) of B r is given in terms of the SVD by 1 

adj B t = £/[diag(0, 0, s { s 2 )]V T . (27) 

We also note that 

BB t B = ^[diag(^^0)]V r . (28) 

These allow us to write 

(A 2 ~ v 2 )fi + ^adj5 r - BB t B = Xs { s 2 UV T = Xs x s 2 A opt , (29) 

where 

A = s x + 5 2 =tra cq(AB t ). (30) 

We can compute the optimal attitude without actually performing the expensive SVD of B if we can find an alternative means 
of computing the quantities appearing in equation (29). Direct computation from equation (17) gives 

adj B t = a,a 2 (b, xb 2 )(r, xr 2 ) r = a,a 2 |b, xb 2 ||r, xr 2 |b 3 r 3 r . (31) 

Then we see from equation (27) that 

5 ,s 2 = ||adj B t \\ f = a x a 2 [b, x b 2 ||r, x r 2 |, (32) 

where ||M|| r denotes the Frobenius (or Euclidean, or Schur, or Hilbert-Schmidt) norm 15,16 

\\M\\ f =[trace(A/Af r )] 1/2 . (33) 

We note from equation (32) that s 2 = 0 if either of the cross products vanishes, as was mentioned above. A little effort is 
required to show that 

A 2 = if +s 2 2 +25,s 2 =||B1|> +2 ai a 2 |b, x b 2 ||r, x r 2 | = a; + a 2 2 +2«,a 2 [(b 1 -b.Xr, r 2 ) + |b, xb 2 ||r, xr 2 |], (34) 

In the two-observation case, A is just the positive square root of the quantity on the right side of equation (34), finding A in 
the case of more than two observations requires solving a quartic equation. To complete the analytic derivation, we need to 
evaluate 

BB t B= £ a,r I; a i b i (i: 'r y )(b J b t )r 1 r (35) 

Combining all these intermediate results with much vector algebra gives the final equation for the optimal attitude estimate: 

\pt s ( fl 1 /A)[b l r, r +(b, xb 3 )(r, xr,) r ]+(a : M)[b 2 r[ +(b,xb J )(r 2 xr,) r ] + b J r, r . (36) 

It is interesting to note that this expression has a unique limit as either or a 2 goes to zero, with A equal to the non-zero 
weight in the limit. This is true even though Wahba’s loss function of equation (15) does not have a unique minimum in 
either limit, since it effectively only includes a single observation. In fact, the limit of the optimal estimate is the TRIAD 
estimate as a~> goes to zero, and A n as a x goes to zero. It is also true, but more difficult to see, that the optimal estimate is 
equal to A n for equal weights, a, = a 2 . 



(37) 


The optimal estimate maps the two reference vectors as 

A v « r . =(«|/ A )^n r i =(a,/A)b, + (a,/A){(r, • r,)b, +[b, - (b, - b, )b, ]|r, x r,|/|b, x b,|} 

and 

=(aJk)A n r, + («:/A)A r2 r 2 =(a,/A)b 2 + (a,/A){( r 1 r,)b 2 +[b, -(b, b 2 )b,]|r, xr,|/|b, xb 2 |}. (38) 

The main point to note about these equations is that the optimal attitude estimate maps both r, and r, into the plane spanned 
by b, and b, . It’s clear from the loss function of equation (15) that this has to be the case; any out-of-plane component 
would be non-optimal. 

In the case that b, b 2 = r, r 2 , equation (34) for A simplifies to A = a, +a 2 , and the optimal estimate is 

A>p< = ( a i^n a 2^T2^/( a i (39) 

Since A n and A n are equal in this case, we see that A opl is equal to their common value, also. 

Mortari has found an alternative representation of the closed-form solution to the two-observation Wahba problem that is 
equivalent to the solution found here l7 . 

OPTIMIZED TRIAD 

Bar-Itzhack and Harman 6 have proposed using equation (37) even when b, -b 2 * r, r 2 . In general, this estimator is not 
optimal, nor is the resulting attitude estimate exactly orthogonal. In order to produce a more nearly orthogonal attitude matrix, 
they employ the first-order orthogonalization step 

A ot =l[(a, + a 2 )” l (a 1 A 7 -| + a 2 A n ) + (a, +a 2 )(a,A^ +a 2 A^ 2 )"‘] (40) 

They call the resulting estimator “Optimized TRIAD.” This estimate has the correct limits of A n and as a, or a 2 tends to 
zero, respectively, but is not the same as A n for equal weights. It avoids the computation of A that is required for the optimal 
estimate, but requires the inverse of a 3x3 matrix. 

There is an alternative way to orthogonalize the matrix computed by equation (37) when b, • b 2 ^ r, • r 2 . This is to extract a 
quaternion from the attitude matrix and then normalize the resulting quaternion. It is well known that the attitude matrix 
computed from a normalized quaternion is guaranteed to be orthogonal 1011,18 . The extraction of the quaternion requires a 
square root, but it is often desirable to compute a quaternion for data transmission or storage, because it stores complete 
attitude information in four components instead of the nine required for the attitude matrix. 

DIRECT QUATERNION METHOD 

All the methods considered so far compute the attitude matrix. If a quaternion is desired, it can be extracted from the attitude 
matrix. However, it would be desirable to avoid this indirect and somewhat costly procedure. Reynolds has proposed a very 
simple estimation algorithm that computes a quaternion directly 7,8 . 

We first present some background information on quaternions to establish our conventions. A more complete discussion can 
be found in reference 1 1 . A quaternion q has a vector part q and a scalar part q s , which we write as 

<7 = [q^J' < 41) 

This is similar to Reynolds’s notation except that we use square brackets rather than parentheses. A unit quaternion (i.e., a 
quaternion with |q| 2 + q] = 1) can be used to represent an attitude matrix, which rotates a vector by 

A(q)\ = {q] - |q| 2 )v + 2(q • v)q - 2q s (q x v) . (42) 

We will follow Shuster’s convention for quaternion products", writing 

p <E> q = [p, p s ] ® [q, q, ] = fop + p s q - p x q, p s q s - P q] ( 43 ) 

This differs from the historical convention in the sign of the cross-product, and has the advantage that the order of quaternion 
multiplication is the same as the order of attitude matrix multiplication: 

A(p® q) = A(p)A(q) . 


(44) 



The quaternion corresponding to the rotation matrix /?(e,0) is 


<1 = 


. 0 0 
esm — ,cos — 
2 2 


(45) 


The derivation of the direct quaternion method begins with the observation that the quaternion that maps the reference vector 
r x into the body frame vector b h using the minimum-angle rotation, is 


<7 mini “ 


V2(I + b, -r,) 


[b, xr,, 1 + b, • r, ] 


(46) 


The most general matrix that maps r, into b { is R(b l ,<t> b )A(q minl )R(r l ,0 r ) * where 0 b and 0 r are arbitrary angles of rotation 
about b! and r u respectively. This general rotation has the quaternion representation 


<7i 


1 


V2( i+M) 

i 


b } sin — , cos — 


®[b, xr,, 1 + bj -rj® 


r, sin — , cos — 
1 2 2 


V 2(l + b,-r 1 ) 


0 0 0 
cos — (bj xr^ + sin — (b 1 + r t ), (1 + bj -r^cos — 


(47) 


where <l>= + <p r . By parallel reasoning, the most general quaternion that maps r 2 into b 2 is given by 


q 2 = rrr cos^(b 2 xr 2 ) + siAb 2 +r 2 ), (l + b 2 -r 2 )cos^ 

-y 2(1 + b 2 • r 2 ) L ^ ^ ^ . 


(48) 


for some angle y/. The vector part of q x is perpendicular to b x - r lf and the vector part of q 2 is perpendicular to b 2 - r 2 . Based 
on this observation, Reynolds proposed to look for a quaternion whose vector part is perpendicular to both b! — and b 2 - r 2 . 
The vector part of q x will be perpendicular to b 2 - r 2 if we choose 

cos(0/2) = ±{[(b, xr,)-(b 2 -r 2 )] 2 + [(b, +r,)-(b 2 - r 2 )] 2 r 1/2 (b, + r,)(b 2 -r 2 ) (49a) 

and 

sin(0/2) = +{[(b, xr,)-(b 2 -r 2 )] 2 +[(b, + r,)-(b 2 -r 2 )] 2 }“ 1/2 (b, xr,)-(b 2 -r 2 ). (49b) 

Substituting this into equation (47) gives 

q x =c,' 1/2 [(b, -r,)x(b 2 -r 2 ), (b, + r,)-(b 2 -r 2 )], (50) 

where c x is the normalization factor 

c, = |(b, - r, ) x (b 2 - r 2 )f + [(b, + r, ) • (b 2 - r 2 )] 2 . (51) 

We have ignored the ambiguous overall sign of the quaternion, which has no significance, since the attitude matrix is a 
homogeneous quadratic function of the quaternion. The appearance of the cross product (bj - r x ) x (b 2 - r 2 ) is not at all 
surprising, since this vector is guaranteed to be orthogonal to both b x - r x and b 2 - r 2 . 

Similarly, choosing yrso that the vector part of q 2 will be perpendicular to bj - r, gives 

q 2 =c 2 ' ,/2 [(b 1 -r,)x(b 2 -r 2 ), (b, +r 2 )-(r, -b,)], (52) 

The vector parts of q x and q 2 are equal up to the normalization constant. However, the scalar part of q x is 

q u =cr ,,2 (b, + r,)-(b 2 - r 2 ) = c' ,/2 [(b 2 -r, -b, r 2 ) + (b, b 2 -r, -r 2 )] (53) 

and the scalar part of q 2 is 

q 2s = c 2 1/2 ( b, +r 2 )-(r, -b,) = c 2 1/2 [(b 2 -r, -b, -r 2 )-(b, -b 2 -r, r 2 )]. (54) 

Thus, q x and q 2 are identical if b, • b 2 = r, • r 2 , and they are equal to 

< 7 3 =c 3 ' 1/2 [(b, -r,)x(b, -r 2 ), b 2 r, - b, r 2 ]. (55) 

We see that q x> q 2y and q 3 all have the same rotation axis, and the rotation angle of q 2 is intermediate between those of q x and 
q 2 . The quaternion q 2 , which treats the two measurements symmetrically, is the estimate preferred by Reynolds; but we will 
also consider the asymmetrical estimates q x and q 2 . 



COMPARISON OF THE DIRECT QUATERNION METHOD WITH TRIAD 

It would seem that the quaternion q x should correspond to the TRIAD estimate A n , q z to /\ n , and q 3 to A n . As evidence lor 
this, we note that the direct quaternion estimation methods have A{q x )r x = A n rj = b t , A(q z )r z = A T2 r z = b 2 , and q 3 
symmetric in the measurements, as A n is. However, we shall now show that this correspondence is not exact. The algebra in 
the general case becomes rather messy, so we consider a simple example. Assume that we have two reference vectors 

r, = [1,0, 0] r and r 2 =[0,l,0f (56) 

and two observation vectors 

b,=[0,0,l] r and b 2 = [cos t?,0,sin t?] r . (57) 

We note that b { • b 2 = r, ■ r 2 only if sin t? = 0, in which case all algorithms should give the same estimate. 


We first compute the TRIAD estimates. Straightforward algebra results in 


0 1 L and A T 


0 10 -sin# cost? 0 -sin(t?/2) cos(t?/2) 0 

A Tl =001, A T2 = 0 0 1 , and A T3 = 0 0 1 . (5 

1 0 0 cost? sint? 0 cos(t?/2) sin(t?/2) 0_ 

We note that A n r { = b t , A r2 r 2 = b 2 , A T3 r + = b + , and A T3 r_ = b_, as expected. However, 

\A n r z - b 2 | = |A r2 r, - bj = 2|sin(t?/2)|, (59 

and 

|^r 3 r i - b .| = \ A Ti r 2 ~ b 2 | = 2|sin(fl/4)|. (59! 

These results are not surprising, since the vectors A n r 2 , A T2 r x , A T3 r v and A T3 r 2 are all in the plane spanned by b t and b 2 , as 
we argued was necessary for an optimal estimator. For comparison with the direct quaternion method, it is interesting to 
present the quaternions extracted from A n , A n , A n : 

q n =i[ 1,1, 1,1], (60 

( 1t2 ~ i[VT -sint?, Vl + sin t?, Vl + sint?, Vl - sint? j , (601 

and 

q n = | [Vl-sin(z?/2), V 1 + sin(i?/2), ^1 + sin(t?/2), ^1 -sin(t?/2)] , (60 

where we have written out the three components of the vector part of each quaternion explicitly. 

The estimates produced by the direct quaternion method embodied in equations (50), (52), and (55) are 

q x = \ (1 + cos t? sin t?)" 1/2 [1, cos t? + sin t?, 1, cos t? + sin t?], (61; 

q 2 = |[l, cost? + sin t?, 1, cost? - sint?], (611 

and 

q 3 = (4 -»- 2 cos t? sin t? - sin 2 t?)" l/2 [l, cost? -(-sint?, 1, cos t?]. (6b 

It is immediately apparent that the quaternions in equation (61) do not correspond to those in equation (60), unless sin t? = 0 
and all reasonable estimators agree. The attitude matrices computed from q [y q ly and q 3 lead to further insights: 




I + cost? sin t? 


0 cost? + sin t? -cost? sint? 

0 cost? sint? cos t? + sint? 

1 + cost? sint? 0 0 

-cost? sint? cost? sin 2 t? 
sint? 0 cost? 

cos 2 t? sint? -cost?sint? 



and 


Mq^) 


1 

4 + 2 cos # sin & - sin 2 # 


-sin #(2eos $ + sin$) 
2sin & 

4 + 2 sin #(cos - sin #) 


2(2 cos & + sin$) 
sin #(2 cos d - sin #) 
2sin & 


-2sin #(costf -sin d) 
2(2 cos # + sin &) 
-sin ^(2 cos & -l- sin #) 


We note that A(q y )r x = b { and A(< 7 2 )r 2 = b 2 , as expected. However, in the general case, 


Vl + cos# sin# \A(q { )r 2 - b 2 | = |A(^ 2 )r, - b t | = |sin #| 
and 

\A(q 3 )r, - b t | = |A(^ 3 )r 2 - b 2 | = V2 (4 + 2 cos # sin # - sin 2 #)" ,/2 |sin #| . 


(62c) 


(63a) 

(63b) 


These residuals are all larger than the corresponding residuals in equation (59), because the vectors A(q x ) r 2 , A(q 2 ) r, f 
/\(< 7 3 )r lf and A(q 2 )r 2 all have components along the y axis in the body frame, which is perpendicular to the plane spanned by 
b x and b 2 . According to our previous argument, they can’t correspond to optimal estimates for any choice of weights. We 
may be prepared to give up optimality for computational simplicity, however. 

SINGULARITY OF THE DIRECT QUATERNION METHOD 

The direct quaternion method has the disadvantage of being ill determined whenever both the vector part and the scalar part 
of the estimated quaternion take the indeterminate value 0/0. We can easily see from equation (50) that q x is undefined if 
b 2 = r 2 , which is when the axis of the attitude rotation is along r 2 (and therefore is also along b 2 ). Similarly, equation (52) 
shows that q 2 is undefined if bj = r lf which is when the axis of the attitude rotation is along r t and b } . These estimators are 
identical in the absence of measurement noise, and we certainly don’t want to depend on measurement noise to avoid a 
singular condition. Thus we see that the direct quaternion method is singular whenever the attitude rotation axis is along r, or 
r 2 (or along or b 2 ). We will now show that the direct quaternion method is singular whenever the attitude rotation axis is in 
the r,, r 2 plane, which means that it is also in the b b b 2 plane. 

If neither b! - r x nor b 2 - r 2 is zero, the vector part of the quaternion estimate vanishes if they are parallel, that is, if 

b 2 -r 2 = fi(b { -rj (64) 


for some scalar /J. The vector b 2 = r 2 + j3(b 1 - Tj) has unit norm, which means that 

l = l+2£r 2 -(b, -r 1 ) + 2/? 2 (l-b 1 -r,). (65) 

Solving this for /? (the zero root is not allowed since b 2 - r 2 * 0) and substituting into equation (64) gives 

b 2 = r 2 -[r 2 -(b, •r 1 )](b 1 -r^. (66) 

It is now straightforward to show that equation (2) is obeyed and that 

b 2 -r 1 =b l -r 2 . (67) 

Thus the vanishing of the vector part of the quaternion estimates of equations (50), (52), and (55) ensures that the scalar parts 
vanish automatically. 

Now let us see what these conditions imply about the attitude quaternion, which certainly exists even if it cannot be 
computed by the direct quaternion method. Equation (42) requires 

b, = {q] - |q| 2 )!■ + 2(q • r, )q - 2q s (qxr,) for/ =1,2. (68) 


From this equation, we can see that 


b 2 r, -b, r, =4^q (r, xr 2 ). 


(69) 


This means that the scalar part of the direct quaternion estimate vanishes either if q is perpendicular to r, x r,, which is to say 
that it is in the r,, r, plane, or else if q , is zero, which indicates a 180° rotation. We still have to investigate the requirement 
that b, - r, is parallel to b, - r,. If q is in the r , , r 2 plane, we can write 


q = a, r, + a, r 3 . 


(70) 



With equation (68), this gives 

b, -r, = 2a,[-(a> + a,r, r,)r, + («, + a,r, r,)r, + < 7 ,(r, xr,)| (71a) 

and 

b : -r, = -2a,[-(a, + a,r, r/)^ +(a, H-c^r, -r,)r 2 +^/ i (r l xr,)]. (71b) 

These two vectors are clearly parallel. On the other hand, equation (68) for a 180° rotation gives 

b, -r, = -2r t +2(q • r,)q = 2q x(q xr,), (72) 

and a straightforward but tedious calculation gives 

(b, - r, ) x (b, -r,) = 4[q-(r, xr,)]q. (73) 


Thus the attitude rotation axis is required to be in the r,. r 2 plane for the 180° rotation case to be singular, also. Thus we have 
completely characterized the singular cases of the direct quaternion method as those cases for which the attitude rotation axis 
is in the r,, r 2 plane, and therefore in the b,, b 2 plane, also. 

The direct quaternion estimate method is singular if the attitude matrix is the identity, giving r, = b, and r 2 = b 2 . We can say 
that the rotation axis is in the r„ r 2 plane in this case, also, because the rotation axis can be arbitrarily assigned for zero 
rotation angle. Reynolds has proposed a method to avoid the singular condition in most cases, but it does not avoid the 
singularity for attitude matrices close to the identity 7 8 . The singular condition can be avoided in all cases by applying 
Shuster’s method of sequential rotations 1019 . This method solves for the attitude with respect to reference coordinate frames 
rotated from the original frame by 180° about the x, y, or z coordinate axis. That is, we solve for the quaternions 

q‘ = <7 ® [e, , 0] = [q, <7 S ] ® [e ■ , 0] = [q s t t -qxe,,-q-e,] for/= 1,2,3, (74) 

where e, is the unit vector along the ( Ul coordinate axis. These rotations are easy to implement on the reference vectors, since 
they simply change the signs of the components perpendicular to e,. Merely permuting and changing signs of the components 
of the rotated quaternion recovers the unrotated quaternion. For example 

q' =[< 7 |. <fc,q 3 ’ < 7 S ]®>[ 1 0,0,0] = [<?,, -93’ (75) 

The method of sequential rotations always avoids the singularity, since the 3x4 matrix 

[q:<? s *i -qxe,:^e 2 -qxe 2 ;q s e 3 -qxe 3 ] (76) 

always has rank three, as can be seen with some effort. Thus the rotation axes produced by the method of sequential rotations 
span the entire three-dimensional space, which means that they cannot all be coplanar with I - ! and r 2 . 

To use Shuster’s rotations to avoid the singularity, we compute the reference vectors r, and r 2 in all four reference frames, 
unrotated and rotated about the x, y, and z axes. We compute the magnitude squared of the cross product (b, - r,) x (b, - r 2 ) 
in each frame, and evaluate the quaternion in the frame where the cross product has the largest magnitude. The above 
analysis shows that this should provide the most robust estimate. If the optimal reference frame is not the unrotated frame, we 
recover the desired quaternion that transforms the unrotated reference frame to the spacecraft body frame by using equation 
(75) or its equivalent for other rotations. 

COMPUTATIONAL EFFORT 

The speed comparison is based on the floating point operation (flop) counts in MATLAB implementations of the algorithms, 
which have the advantage of being platform-independent. There are some caveats to make with regard to timing comparisons. 
First, for ground computations, absolute speed isn’t all that important, since the estimation algorithm is only a part of the 
overall attitude determination data processing effort. Speed was more important in the past, when thousands of attitude 
solutions had to be computed by slower machines. Second, for real-time processing, as for an attitude control system onboard 
a spacecraft, the longest time is more important than the average time, because the attitude control system processor has to 
finish its task in a limited amount of time. This works against methods that may require sequential rotations. 

Four methods for computing the attitude matrix are compared in Table 1 : asymmetric TRIAD of equation (8), symmetric 
TRIAD of equation (12), the optimal two-measurement estimator of equation (36), and Optimized TRIAD of equation (40). 

An “asymmetric” estimator maps one of the two reference vectors into the corresponding observed vector exactly, throwing 
all the measurement errors into the other vector. A “symmetric” estimator, on the other hand, treats the two measurements 
symmetrically. The cost of using these four estimators to produce a quaternion is also presented. Every algorithm except 



Optimized TRIAD computes the quaternion by extracting it from the corresponding attitude matrix, a process that costs 29 
flops (see the Appendix). The quaternion output of Optimized TRIAD is cheaper than the attitude matrix output because it is 
extracted from the estimate of equation (39) rather than trom equation (40). In addition to these tour estimators, three other 
estimators are included for quaternion output only: the asymmetric direct quaternion estimator of equation (50), the 
symmetric direct quaternion estimator of equation (55), and QUEST, for comparison The computational effort for the 
direct quaternion estimation methods is given both with and without the use of rotations to avoid singular configurations. The 
computational effort of QUEST does not include the cost of sequential rotations. No special efforts have been made to 
achieve the most efficient possible implementation of any of the algorithms. 


Table 1: Computational Effort of Estimation Algorithms in Flops 


Algorithm 

A output 

q output 

Asymmetric TRIAD 

143 

172 

Symmetric TRIAD 

166 

195 

Optimal Two-Measurement Estimator 

265 

294 

Optimized TRIAD 

335 

273 

Asymmetric Direct Quaternion 


46 

Asymmetric Direct Quaternion with Singularity Avoidance 


108 

Symmetric Direct Quaternion 


50 

Symmetric Direct Quaternion with Singularity Avoidance 


112 

QUEST 


190 


Several conclusions are apparent from these results. Symmetric estimators are a little more expensive than asymmetric 
estimators, in general. Optimized TRIAD with the approximate matrix orthogonalization of equation (40) is significantly 
more expensive than the optimal two-measurement estimator. If quaternion output is desired, Optimized TRIAD is slightly 
less expensive than the optimal two-measurement estimator; but the savings are less than 10%. However, the optimal two- 
measurement estimator and Optimized TRIAD (and even symmetric TRIAD) require more computational effort than QUEST 
to produce a quaternion. Asymmetric TRIAD is only slightly less expensive than QUEST, but the direct quaternion 
estimation methods developed by Reynolds are significantly faster. The implementation of rotations to avoid singularities in 
the direct quaternion estimation methods more than doubles their computational cost, but they are faster than other methods 
even with this modification. None of the three algorithms faster than QUEST is optimal, though; and QUEST also has the 
advantage of being a general-purpose algorithm applicable to any number of measurements, which avoids the need to 
develop and test a special-purpose two-observation algorithm. 


ACCURACY 

We will analyze two test scenarios, using the nine estimators with quaternion output that were used in the timing tests. The 
first scenario simulates two star trackers with narrow fields of view and orthogonal boresights at [1, 0, 0] and [0, 1,0]. We 
shall assume that the first tracker is tracking five stars at 
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We simulate 1000 test cases with random attitude matrices. We use the attitude matrices to map the eight observation vectors 
to the reference frame, add Gaussian random noise with equal standard deviations of 6 arcseconds per axis to the reference 
vectors, and then normalize them. The errors are unconventionally applied to the reference vectors rather than the observation 
vectors so that equation (77) will remain valid in the presence of noise. The two-observation estimators use averages of the 
multiple vectors observed by each tracker, as suggested by Bronzenac and Bender 5 . In this example the two averaged vectors 
in the body frame are along the star tracker boresights. The optimal estimator weights for these estimators are proportional to 
the inverse measurement variances, or to the number of vectors included in the average, so we use a 2 = 0.6 a { tor the optimal 
two-measurement estimator and Optimized TRIAD. 

We treat the eight star measurements independently in QUEST, rather than averaging them. QUEST requires 316 flops for 
eight measurements, but avoids the cost of averaging the vectors, which is 108 flops. Thus the computational effort of 
QUEST should be taken as 208 flops for comparison with the other estimators in this eight-measurement example, making it 
more expensive than the direct quaternion estimator and TRIAD, but faster than the optimal two- vector estimator and 
Optimized TRIAD. In these tests, QUEST always used information about the true quaternion to determine the optimally 
rotated reference frame for estimation, without the need to perform sequential rotations. 

Table 2 shows that symmetric TRIAD, the optimal two-measurement estimator, and Optimized TRIAD perform almost as 
well as QUEST. This justifies Bronzenac and Bender’s procedure of using average observation and measurement vectors for 
the two star trackers. It should be noted, however, that the choice of orthogonal tracker boresights is optimal for this 
approximation, and that symmetric TRIAD is the only one of these algorithms that is computationally less expensive than 
QUEST, requiring 13 fewer flops. The symmetric direct quaternion estimator with singularity avoidance provides average 
and maximum errors within 10% of those of the best estimators with less computational effort, though. 


Table 2: Average (Maximum) Estimation Errors (arcseconds) for Star Tracker Scenario 


Algorithm 

All Cases 

w 

IV 

To 

1^1 <1/2 

Asymmetric TRIAD 

4.6(12.1) 

4.5 (11.3) 

4.7 (12.1) 

Asymmetric Direct Quaternion 

13.6(2562) 

5.2(17.8) 

20.1 (2562) 

Asymmetric Direct Quaternion with Singularity Avoidance 

5.1 (16.9) 

5.1 (14.5) 

5.1 (16.9) 

Symmetric TRIAD 

4.4(12.2) 

4.3 (11.6) 

4.5 (12.2) 

Symmetric Direct Quaternion 

14.2 (4763) 

4.7 (14.6) 

21.6(4763) 

Symmetric Direct Quaternion with Singularity Avoidance 

4.7(12.9) 

4.6(12.1) 

4.8 (12.9) 

Optimized TRIAD 

4.6(12.1) 

4.5 (11.3) 

4.7(12.1) 

Optimal Two-Measurement Estimator 

4.6(12.1) 

4.5 (11.3) 

4.7(12.1) 

QUEST 

4.4(11.8) 

4.3 (11.5) 

4.4 (11.8) 


The results also show that symmetric estimators perform slightly better than asymmetric estimators in this scenario. This was 
expected, since the number of stars tracked in the two trackers and thus the measurement weights are nearly equal. A 
symmetric estimator would be preferred in a real star tracker application, since there would be no way of predicting a priori 
which tracker would view more stars. 

Table 2 also shows inferior performance of the direct quaternion estimators without singularity avoidance. The performance 
is not so bad in the 436 simulated cases with \q 3 \ > 1/2 as in the 564 cases with \q 3 \ < 1/2 . The latter are the cases in which we 
would expect singularities to occur, since they have either small rotation angles or rotation axes close to the x-y plane, the b |t 
b 2 plane in this scenario. This shows the importance of avoiding singular cases in an application of these estimators. We note 
that the performance with singularity avoidance, as well as the performance of all the other estimators, is independent of q 3 . 

The second scenario that we consider is a sun-mag system, similar to that on SAMPEX 3 , assuming a digital sun sensor with 
accuracy of 0.1° and a magnetometer with effective accuracy of 1°. We assume that the Sun is at the center of view of the 
digital sun sensor at bj = [ 1, 0, 0] r , but the orientation of the magnetic field vector is not fixed in the spacecraft body frame. 
We simulate 1000 random attitude matrices and random magnetic field vector orientations, except that we do not allow the 




magnetic field direction to be within 5° of the ± v axis. These are the cases with coaligned vectors that the SAMPEX onboard 
attitude determination system rejects. We use the attitude matrices to map the Sun and magnetic field observation vectors to 
the reference frame, add Gaussian random noise w'ith standard deviations as specified above, and then normalize the 
reference vectors. In this case the optimal estimator weights have a 2 = 0.0 1 

The estimation errors for this scenario are presented in Table 3. The roll (x axis) and pitch/yaw (root-sum-squared of v and z 
axes) errors are presented separately, since the estimate of pitch and yaw provided by the digital sun sensor on the x axis is 
more precise than the roll angle estimate provided by the magnetometer. We note from these tables that QUEST and the 
optimal two-measurement estimator give identical errors, as they must since this scenario has two vector measurements. 
Since the weight assigned to the magnetometer measurement is so much less than the weight of the sun sensor measurement, 
Optimized TRIAD and asymmetric TRIAD give virtually the same results as the optimal estimators. The asymmetric direct 
quaternion estimator with singularity avoidance provides equivalent pitch and yaw errors, and average and maximum roll 
errors within 20% of those of the best estimators, with less computational effort. 

Symmetric estimators are inferior to asymmetric estimators in the sun-mag scenario, since they allow the magnetometer 
errors to corrupt the sun sensor determination of pitch and yaw. Table 3a shows that the direct quaternion estimation method 
must be modified to provide acceptable roll estimation in the 551 cases with |qj < 1/2, where qjs the component of q 
perpendicular to the bj, b 2 plane. Table 3b shows that pitch and yaw estimates provided by the asymmetric direct quaternion 
estimator are insensitive to these singular configurations, since this estimator maps r t into bj exactly. 


Table 3a: Average (Maximum) Roll Estimation Errors (degrees) for Sun-Mag Test Case 


s : : 

Algorithm 

All Cases 

M^/2 

|qJ<i/2 

Asymmetric TRIAD 

0.88 (3.06) 

0.93 (2.92) 

0.84 (3.06) 

Asymmetric Direct Quaternion 

2.82(114) 

1.07 (4.78) 

4.24(114) 

Asymmetric Direct Quaternion with Singularity Avoidance 

1.01 (3.58) 

1.05 (3.36) 

0.98 (3.58) 

Symmetric TRIAD 

0.88 (3.06) 

0.93 (2.92) 

0.84 (3.06) 

Symmetric Direct Quaternion 

1.97 (86.5) 

0.98 (3.78) 

2.77 (86.5) 

Symmetric Direct Quaternion with Singularity Avoidance 

0.92 (3.23) 

0.96(3.16) 

0.89 (3.23) 

Optimized TRIAD 

0.88 (3.06) 

0.93 (2.92) 

0.84 (3.06) 

Optimal Two-Measurement Estimator 

0.88 (3.06) 

0.93 (2.92) 

0.84 (3.06) 

QUEST 

0.88 (3.06) 

0.93 (2.92) 

0.84 (3.06) 


Table 3b: Average (Maximum) Pitch/Yaw Estimation Errors (degrees) for Sun-Mag Test Case 


s i : 

Algorithm 

All Cases 

Nxl^l/2 

kl<i/2 

Asymmetric TRIAD 

0.13 (0.36) 

0.13(0.35) 

0.13(0.36) 

Asymmetric Direct Quaternion 

0.13 (0.36) 

0.13 (0.35) 

0.13 (0.36) 

Asymmetric Direct Quaternion with Singularity Avoidance 

0.13(0.36) 

0.13 (0.35) 

0.13(0.36) 

Symmetric TRIAD 

0.43 (1.60) 

0.42 (1.54) 

0.43 (1.60) 

Symmetric Direct Quaternion 

1.53 (96.3) 

0.50(1.91) 

2.37 (96.3) 

Symmetric Direct Quaternion with Singularity Avoidance 

0.48 (1.92) 

0.49 (1.91) 

0.48 (1.92) 

Optimized TRIAD 

0.13 (0.37) 

0.13 (0.35) 

0.13 (0.37) 

Optimal Two-Measurement Estimator 

0.13 (0.37) 

0.13 (0.35) 

0.13(0.37) 

QUEST 

0.13 (0.37) 

0.13 (0.35) 

0.13 (0.37) 






























CONCLUSIONS 


We have analyzed four spacecraft attitude determination methods using exactly two vector measurements: the well-known 
TRIAD algorithm, an optimal closed-form two-measurement of Wahba’s optimization problem, the Optimized TRIAD 
algorithm of Bar-Itzhack and Harman, and the direct quaternion estimation method of Reynolds. These methods are 
applicable to a variety of problems, including coarse “sun-mag” attitude estimation using the unit vector to the Sun and the 
Earth's magnetic field vector and precise estimation using unit vectors to stars tracked by two star trackers. For TRIAD and 
the direct quaternion estimation method, we investigate both “asymmetric” forms that map one of the two reference vectors 
into the corresponding observed vector exactly, and “symmetric” forms that distribute the errors symmetrically between the 
two measurements. We also include the well-known QUEST algorithm for comparison, 

The computational speed of the algorithms was compared using floating point operation (flop) counts in MATLAB. These 
show that Optimized TRIAD and the optimal two- measurement estimator are more expensive than QUEST, which has the 
additional advantage of being a general-purpose algorithm applicable to any number of measurements. The direct quaternion 
estimation methods are significantly faster than QUEST, however. Both QUEST and the direct quaternion estimation 
methods have the disadvantage of sometimes requiring special computations to avoid singular cases, but the direct quaternion 
estimation methods are faster than any other methods even with these modifications. 

We analyzed the accuracy of the estimators in two test scenarios. The first scenario simulated two star trackers with narrow 
fields of view and orthogonal boresights, using average vectors for five stars in the first tracker and three in the second. The 
second scenario simulated a digital sun sensor with accuracy of 0.1° and a magnetometer with effective accuracy of 1°. 
Symmetric estimators outperformed asymmetric estimators in the first scenario, and asymmetric estimators were superior in 
the second, as was expected. With this proviso, all the estimators had comparable errors. The one exception is that the direct 
quaternion estimators had larger errors if not modified to avoid singularities, showing the need for these modifications. 

This paper demonstrates the superiority of TRIAD, QUEST, and the direct quaternion estimation methods for attitude 
estimation from two vector measurements. 
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APPENDIX 

The following MATLAB function was used to extract quaternions from attitude matrices 1011 , 
function q = dcm2quat(a) 

% finds the quaternion representation of a direction cosine matrix a 

% find maximum of trace or diagonal element of direction cosine matrix 
tra = trace (a) ; 

[mx,i] = max ( [a ( 1 , 1 ) a (2, 2) a (3, 3) tra]); 

% compute unnormalized quaternion 

if i==l, q = [ 2*mx + l-tra; a ( 1 , 2 ) +a (2 , 1 ) ; a ( 1 , 3) +a ( 3 , 1 ) ; a { 2 , 3) -a ( 3 , 2 ) ] ; end; 

if i“2, q = [a (2, 1) +a (1, 2) ; 2*mx+l-tra; a (2, 3) +a (3, 2) ; a { 3 , 1 ) -a ( 1, 3 ) ] ; end; 

if i==3 , q = [a(3,l)+a{l,3);a(3,2)+a(2,3) ; 2*mx+l-tra; a (1, 2) -a (2, 1) ]; end; 

if i— 4, q = [a(2, 3)-a{3,2) ;a(3, l)-a(l,3) ;a(l,2)-a(2, 1) ;l+tra] ; end; 

% normalize the quaternion 
q = q/norm (q) ; 



